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Abstract 



This thesis covers various aspects of open systems in classical and quantum mechanics. In the first 
part, we deal with classical systems. The bath-of-oscillators formalism is used to describe an open 
system, and the phenomenological Langevin equation is recovered. The Fokker-Planck equation is 
derived from its corresponding Langevin equation. The Fokker-Planck equation for a particle in 
a periodic potential in the high-friction limit is solved using the continued-fraction method. The 
equilibrium and time-dependent solutions are obtained. Under strong periodic driving, we observe 
significant non-linear effects in the dynamical hysteresis loops. Shapiro steps appear in the time- 
average of the drift velocity curves. Similar study is carried out for a dipole in an electric field. 

In the second part of the thesis, we begin the study of open quantum systems by re-deriving the 
quantum master equation using perturbation theory. The master equation is then applied to the 
bath-of-oscillators model. The subtleties and approximations of the master equation are discussed. 
We then use the master equation to solve the damped harmonic oscillator. The equilibrium solution 
coincides with the canonical distribution. The steady state response to DC and AC forces is also 
studied. Driven systems are more challenging in quantum open systems, and we manage to solve the 
quantum master equation with the continued-fraction method. We obtain the frequency-dependent 
susceptibility curves, which exhibit typical absorption and dispersion profiles. 
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Chapter 1 

Introduction 



Most coffee lovers have the frustrating experience of having their coffees cooled down before they 
could finish them. High school physics tells us that it is due to the heat transferred to the surround- 
ings. In the language of statistical mechanics, it is because of the interaction with the environment. 
Every object, big or small, classical or quantum, is subject to this interaction. 

The focus of this thesis is to study the effects of the environment on the statics and dynamics 
of our systems. The effects are two-fold: fluctuation and dissipation. Think of the pollen grains 
in water as observed by Robert Brown. Their trajectories exhibit random behavior due to the 
collisions with the water molecules. This randomness makes us unable to make exact predictions of 
the system evolution, we can only talk about its statistical properties. During the collisions, some 
of the pollen grains' momentum is transferred to the medium, causing them to lose energy. Due of 
this dissipative process, the system can relax to a stationary state. These random and dissipative 
effects apply to any open systems, and we will address them in both classical and quantum regimes. 

This thesis is organized as follows. We deal with classical open systems in Chapter 2 to Chapter 
4, while Chapter 5 and Chapter 6 are devoted to the study of open quantum systems. In Chapter 
2, we start from a microscopic viewpoint of classical open systems, and introduce two equivalent 
approaches to study them: the Langevin equation (trajectory) and the Fokker-Planck equation 
(distribution). We then solve the Fokker-Planck equation of a particle in periodic potential, and 
study its steady-state properties in time independent and time-dependent fields (Chapter 3). Sim- 
ilar study is done for a classical dipole in Chapter 4. 

In Chapter 5, we start the exploration of quantum open systems by discussing the reduced de- 
scription of the open systems. Then using perturbation theory, we present a concise derivation of 
the so-called master equation: the equation of motion of the reduced density matrix. In Chapter 6, 
we use the master equation to study a damped quantum harmonic oscillator and make comparison 
with exact results when available. 
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Chapter 2 

Classical Open Systems: Langevin 
and Fokker- Planck Equations 



2.1 Introduction 

In 1908, the French physicist Paul Langevin proposed a modified version of the Newton equation 
to describe the dissipative and random behavior of a Brownian particle pp. Consider a particle in 
one dimension for notational simplicity, the phenomenological equation of motion reads 

Mx + ~jMx + V'{x) =£{t). (2.1) 

The term "fMx causes dissipation, and 7 is the damping rate. The random force £(t) originates from 
the impacts with the fluid molecules. It is assumed that the random force is Gaussian distributed, 
and thus fully described by the first two moments: 

<£(*)> = 0, (£(*)£(*')> = 2M 7 k B T5(t - 1'). (2.2) 



Equation (2.1) has a time-local frictional term, thus assuming the friction does not depend on 
the velocity of the past. In many cases of interest, the bath has a finite memory time, and the 
coresponding equation of motion is 

Mx + M f ~t(t-t')x(t') dt' + V'(x) =£(*), (2.3) 
«/ — 00 

which is called the generalized Langevin equation. The damping rate is replaced by a memory 
function j(t), which comes from the finite noise correlation time = 2MksT'j(t — if). 



There are problems with these phenomenological equations: they are not time reversal-invariant 
and therefore contradict Newton's reversible laws. Furthermore, we cannot carry forward this 
formalism into the quantum regime since we do not know how to quantize a dissipative equation. 
In the next section, we will start from a "microscopic" point of view to describe fluctuations and 
dissipation using a time- honored framework in classical physics — Hamiltonian mechanics. 



2.2 Bath-of-Oscillators Formalism 

Here we will derive the Langevin equation by considering the Hamiltonian equations of the global 
system, system plus bath (from now on we will call the system of interest the "system", the sur- 
roundings the "bath"). 
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To obtain the equations of motion, we need the Hamiltonian 



fltot — H sys + //bath + Vint- (2-4) 

This total Hamiltonian //tot consists of the system Hamiltonian H sys , the bath Hamiltonian //bath 
and the interaction term Vi n t • We need to give content to the bath and coupling Hamiltonians. 

Following Zwanzig [2] and others [31 HIE], the bath is modeled as a set of harmonic oscillators. 
The replacement of the true bath by a set of harmonic oscillators is effectively equivalent to the 
assumption that the coupling is weak such that the bath is only slightly perturbed away from its 
equilibrium configuration [3J. By the same argument, it is reasonable to assume that the system- 
bath coupling is linear with respect to the bath coordinates. We write the Hamiltonians as 

2 

#s ys = T^ + Vix), (2.5) 



N 



//bath = E(^- + ^ m ^')' ( 2 - 6 ) 
ot=l a 

N 

V^t = -F{x,p)^c a x Q + AV(x,p), (2.7) 

a=l 

where a denotes the bath modes (which can be continuous), c a the coupling constants. F(x,p) can 
be any function of the system's coordinate and momentum (x,p). It is worth noting that though 
the influence of the system on the bath is small, the opposite might not be true. 

2.2.1 Counter- Term 

An additional term AV is added to the coupling Hamiltonian to compensate the re-normalization 
caused by the term Fx a [6]. In other words, we want to ensure that the global minimum of the 
Hamiltonian is determined by the bare potential V{x) alone. To look for the minimum with respect 
to the bath, we need 

d^ tot = m a u 2 a x a - c a F = 0, (2.8) 

and obtain x a = m c< ^2 F. Using this result, we find the minimum with respect to the system 
coordinate 

&ff tot dV ^ cl dF „ , dAV 



dx dx ^—i raaUj"^ dx dx 



a=l 



In order to satisfy 8 ffi ot = §p, we need the counter-term 



dx dx 



N 2 

^) = E ^2 (2-10) 



a=l 
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2.2.2 Bath-of-Oscillators Hamiltonian 



Finally gathering all the above results, we can write the total Hamiltonian in the form 



tot 



-f^sys + -f^bath + Vint 



(2.11) 



H s y s + 



a=l 



2m f 



This model is frequently called the Caldeira-Leggett model in the context of quantum dissipative 
systems [31 Hj. This framework can also be used to study problems involving spin, in which the 
coupling is a function of spin variable, F = F(s). In the next section, we will show that the 



Hamiltonian (2.11) describes dissipation and fluctuations in the system. 



2.3 The Equation of Motion (Langevin) 
2.3.1 Derivation 

Consider a dynamical variable A(x,p) that only depends on the system's coordinate and momentum, 
the Hamiltonian equations of A(x,p), x a and p a are given by 



dA 
~dt 



N n 
C 

1 V 



N 



a=l 



dt m, 
where the Poisson bracket is 



{A, H sys } + £ — ^ {A, F 2 } - £ c a x a {A, F}, 

dp a 

dt 



{A,B} 



dA dB dA dB 



dx dp dp dx 

The solution to the bath coordinate is that of a forced harmonic oscillator 

-t 



x a (t) = x h a {t) + 



ds sin[w a (t — s)]F(s), 



where the homogeneous part is the free harmonic oscillator evolution 

Xa(t) = x a (t ) cos[u a (t - t )] + Pa ^ sm[u a (t - t )}. 



m a uj a 



(2.12) 



(2.13) 



(2.14) 



(2.15) 



(2.16) 



Substituting Eq. (2.15) and Eq. (2.16) into Eq. (2.12) and performing integration by parts, we 
obtain 



^ = { A, H sys } — {A,F} [^(t) — M ds j(t — s) ^ 



where 



N 



^— ' L m a iof 

a=l 



°^F(t ) cos[w Q (i - t )] 



7(*) 



1 N 

-Y 



M ^— ' m„uj 



cos(u a t). 



(2.17) 

(2.18) 
(2.19) 



a=l 
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We have obtained an equation similar to the Langevin equation for a general dynamic variable 



A(x,p). The first term in Eq. (2.17) is the free evolution. The second term gives rise to the 



fluctuations. The integral term keeps the memory of the previous states and causes dissipation. We 
will justify the interpretation of fluctuations and dissipation in the next subsection. 



2.3.2 Fluctuation and Dissipation 
Fluctuation 

The term £(t) contains the free evolution of the bath oscillators and the initial state of the system. 
Being the sum of many terms with different frequencies and phases, it behaves like a random force 



(see Figure 2.1). Let us assume the initial distribution of the bath follows the classical canonical 
distribution 



JV 



Pb 



ath(*o) = Z 1 exp | - p ^2 [ Pa r!i°^ + \rn a iJ 1 a [x n (/„ i 



a=l 



2m n 



-F(t c 



}■ 



(2.20) 



Drawing x a (to) and p a (to) from the canonical distribution, becomes a random force with 
Gaussian distribution. Then it is fully characterized by the first two moments 



= o, mm) = K(t-t'), 

where the correlator is related to the damping kernel as K(t) = 2M ksTjit). 



(2.21) 



Dissipation 

If we particularize the equation of motion to the system Hamiltonian, the free evolution is zero and 
we are left with 



dH. 



sys 



dt 



-{H ays ,F}Z(t) + {H sys ,F} ! ds 7 (t-a 

Jt 



dF(s) 
ds 



(2.22) 




Figure 2.1: Computer simulation of the noise term, It is obtained from summing (2.18) over 



1000 oscillators having the canonical distribution as initial conditions. It does indeed look random. 
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The first term averages to zero, loosely speaking. Otherwise, we assume T = 0, where there is no 
thermal fluctuation and the first term is automatically zero. We are then left with the integral. For 
simplicity, we consider a free particle H sys = p 2 /2M with coordinate coupling F = x, and delta 
damping kernel j(t) = j5(t), one finds 

— oc - 7P . (2.23) 



The negative rate of change indicates that dissipation indeed takes place. Eq. (2.21) relates the 



noise or fluctuating force to the dissipation, and is known as the fluctuation-dissipation theorem. 
Spectral Density 

It is convenient to introduce the spectral density 

N 2 

J{u) = lY J ^^5(ui-Lj a ). (2.24) 

a=l 



The correlator (2.21) and the damping kernel can then be written as 



Kit) = 2k B T f^— ^M C os(wt); (2.25) 
Jo 71 ^ 

7 * = T7 / ^cosM). (2.26) 

M J IT Ld 

For a bath with discrete modes, the spectral density is a set of delta peaks. However in a dissipative 
bath, the eigenfrequencies oj a form a continuum and J{oj) becomes a smooth function of u. 

All the effects of the bath are incorporated into J(uj), which involves the frequencies and cou- 
plings. It these are not fully known, one proceeds to "model" the bath, assuming different functional 
dependences. In the next section, we will discuss a system (Rubin model) where J(oj) can be ex- 
plicitly computed, and this will provide us insights how to properly do such modeling. 

2.3.3 Examples: Brownian Particle and Spin 



In this section, we study the Langevin equations (2.17) for a Brownian particle (translational motion) 
and a Brownian spin (rotational motion). 

Brownian Particle 

Consider a particle with Hamiltonian H sys = p 2 /2M + V(x), with its coordinate (F = x) coupled 



to the bath. Equation (2.17) for A = x and A = p then gives 



^ = P/M; f t = -V'{x) + ^(t) ~ f ds 7 (t - s)p(s), (2.27) 

which is exactly the generalized Langevin equation introduced phenomenologically at the beginning 
of this chapter. 
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If the bath is Ohmic (Markovian limit), namely J(co) = Mjuj, the damping kernel becomes a 
delta function and we recover the usual Langevin equation 



dp 
dt 



-V'(x) + ~ IP, 



with the bath correlator 



WW)) = 2M 1 k B T5{t-t'). 



(2.28) 



(2.29) 



Brownian Spin 

Equation (2.17) is also valid for ci spin with Hamiltonian H S y S — H S ys[s X i Sy^ 

Using the spectral 

density J (to) = Xu, we write directly the memoryless equation of motion for s 



j t = {s,# sys } + {s,F}[£(t) + A^ 



(2.30) 



s x B c g + s x £(t) — As x A 



ds 
dt'' 



where 



B e flf 

A 



dH< 



sys 



ds 



dF 

~ds~ 
/dF\UdF 
V ds ) V ds 



(2.31) 
(2.32) 
(2.33) 



This is the Lagenvin equation for a spin. The first term arises from the free Hamiltonian and 
causes the spin to precess around the field direction. For example, the Hamiltonian of a spin s in a 
magnetic field B is H sys = — s-B and B e g = B. The second and third terms describe the fluctuation 
and dissipation as discussed before. The damping term is a generalization of the phenomenological 
equations proposed by Gilbert and Landau-Lifshitz [7j. 



2.3.4 Rubin Model |H [S] 

Here we look at an instructive example where the oscillator-bath model represents the actual Hamil- 
tonian. In the Rubin model (see Figure 2.2), a heavy particle of mass M and coordinate x is bi- 
linearly coupled to a half infinite chain of harmonic oscillators with mass m and spring constant 



/ 



r 2 -. The total Hamiltonian is 



tot 



2 00 2 

K +V(x) + J2\^: + ^n+i 



2M 



n=l 



12m 



+ 2 {X 



x ij 



(2.34) 



To cast the above Hamiltonian into the standard form (2.11), we make the following transfor- 
mation to normal modes X(k) 



x n = \/2/ir / dksin(nk)X(k). 
Jo 



(2.35) 



The underlaying canonical variables are x = <p and p = s z 7 . But our formalism does not depend on the 
canonical variables, as we derive the equation of motion for any A(x,p) and we can set A = Sj for i = x, y, z. 
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^vw\aO^ • • • 

M f m 

Figure 2.2: Pictorial Sketch of the Rubin model, figure taken from [5]. 

In normal mode representation, the Hamiltonian reads 

fltot = 7^7 + V(x) + f -x 2 + ^ dk ( P ^ + ^u 2 (k)X 2 (k) - x c(k)X(k)) . (2.36) 
2M 2 Jq \ 2m 2 / 

The eigenfrequency oj{k) and the coupling function c(k) are 

u(k)=oj R \sin(k/2)\; c(k) = J- sm(k). (2.37) 

V 7T 4 



Comparing with Eq. (2.24), the spectral density is found to be 



7T r c 2 (k) muj R ( w 2 \V2 



J(w) = -/ dfc— ^-5[w-w(fc)] = — ^ fX 2- 1 Q(u, R -u), (2.38) 

2 J muj{k) 2 \ uj r / 

where Q(uj) is the Heaviside step function. The frequency ojr is the highest frequency in the bath 
and cuts off the spectral density J(ui) (see Figure 2.3). In fact, in any physical system, there always 
exists a cut-off frequency such that the contribution at high frequencies is suppressed. With the 



above spectral density, the damping kernel Eq. (2.25) becomes 



l(t) = 2M^«— T— . ( 2 -39) 

where J\{uj R t) is the first order Bessel function. The memory time in the kernel is of the order of 
1/ojr. For large ojr (corresponds to a stiff spring), the kernel becomes a sharply peaked function 
around zero (see Figure 2.4). 

The Rubin model not only gives us a concrete example where the oscillators represent a true 
bath, it provides us useful insights about the properties that are common to any physical bath. 
The continuum limit of the spectral density comes from the infinite number of degrees of freedom 
of the chain of oscillators. There is a natural cut-off to the spectral density at high frequency. We 
also see that, at large ojr, we approach the Ohmic limit [Markovian, J(u) oc uj], the damping kernel 
becomes short-lived and the harmonic chain exhibits little retardation. 



2.4 Fokker-Planck Equations 

We have studied the trajectory approach to the open systems based on the Langevin equation. Here 
we will look at the phase space distribution function and its evolution equation, the Fokker-Planck 
equation. Both approaches are equivalent, provided the noise is delta correlated and has a Gaussian 
distribution; we will see why soon. 
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Spectral Density of the Rubin Model 



0.25 




■ 0.15 



0.05 



(D/av 



Figure 2.3: The spectral density of the Rubin model, mw^ = 1. 



The Damping Kernel of Rubin Model. 



0.25 



0.15 




0.05- 



-0.05 



10 12 14 

(0 R t 



20 



Figure 2.4: The damping kernel of the Rubin model, m = M = 1. If ujr is large, the curve is 
sharply peaked. 



2.4.1 Derivation a la Zwanzig [2] 

To derive the Fokker-Planck equation, we start with a general Langevin equation of a set of variables, 
a = {a?}, the equation of motion in vector form is 

da 

- = v(a) + *(*), (2.40) 
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where the noise term £(t) is Gaussian distributed and has the following properties 

<*(*)> = 0, «(t)e(0 T >=2Bd(t-0- (2-41) 

We are interested in the probability distribution of the dynamical variables, w(a, t). The nor- 
malization condition requires 

(2.42) 



(2.43) 



daw(a, t) = 1, for all t. 
Similar to the conservation law in electromagnetisrrj^J we have the continuity equation 



dw d 
+ 



( 



da 



dt da \ dt 



w 



0. 



Replacing ^ by Eq. (2.40), one has 

dw 
~dt 



Cw 



d_ 
da 



v(a)w + £(t)w 



(2.44) 
(2.45) 



It is still a stochastic differential equation since it contains the noise. We are interested in the noise 
average of w(a,t). We denote the noiseless part of the operator 

A><& = • (v(a)$. 



da 



C w 



d 



so that 

dw 

dt da 

The formal solution to the differential equation is (after setting initial time to = 0) 



w(a, t) 



JC 



w(a, 0) 



o 



t Q 

dse^ Co ^--^(s)w(a,s). 
da 



(2.46) 



(2.47) 



(2.48) 



Note that w(a, t) only depends on the noise £(s) at earlier time s < t. Substituting the above 

d 



equation into Eq. (2.47), we obtain 

dw 
~dt 



C w(a,0) - £(t)e tCo w(a,0) 



(2.49) 



d_ 
da 







t a 

dse (t-s)c » .^( s ) w ( a)S ). 
da 



Now we take the average over noise. The second term averages to zero. Since the noise is 
Gaussian, we can express any moments in terms of the first two moments. The last term contains two 
explicit noise factors, £(s) and £(t), and also the implicit noise factors in w(a, s). Since w(a, s) only 
depends on the noise at earlier time s' < s ; the pairing with either £(s) or £(t) gives zero contribution 
as the noise is delta correlated. We only need to consider the pairing (£(i)£(s)) = 2B<5(£ — s). 
Denoting the noise average of the distribution as (w(a,t)) = W(a,t), the result is the Fokker- 
Planck equation 

dW(a,t) d 
dt da 

The first term on the right hand side is what one has in the absence of noise. The effect of the 
noise is introduced by the second term. This term has a diffusion structure, with a second order 
derivative, as in the standard diffusion equations. 

2 In electromagnetism, t 
charge density respectively. 



v(a)W(a,t) 



d - d 



(2.50) 



2 In electromagnetism, we have the continuity equation: V • J + §f 



0, where J and p are the current density and 
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2.4.2 Applications to Brownian Particle and Spin 

As in the previous section, we will see the examples of a Brownian particle and a Brownian spin. 



Brownian Particle 

According to the recipe above, we can write down the Fokker-Planck equation for a Brownian 
particle from its Langevin equation (2.28). The quantities that enter into the general Fokker-Planck 
equation (2.50) are 







v(a) 



p/M 
-V'(x) — jp 



B 





M-/k B T 



The resulting Fokker-Planck equation is called the Klein-Kramers equation [9] 



dW 



dt 



Mdx v ' dp 



d / 

7— [p + Mk B T- 



W. 



(2.51) 



d_ 

dp \ 1J D dp/ 

The first two terms arise from the Liouville operator {H sys , W}, while the last two terms capture 
the effects of the interaction with the bath: damping and diffusion. 



Brownian Spin 

The corresponding Fokker-Planck equation for a Brownian spin is [7] 



dW(s,t) d 



{s x B cfr - As x A [s x (B cff - fc B T— )] }\V(s, t). (2.52) 



dt ds 

We will return to this type of orientational diffusion equation in Chapter 4 when we study the Debye 
dipole. 

2.4.3 Solving the Fokker-Planck Equations [9] 

In most cases, the Fokker-Planck equation is not solvable analytically, and we have to resort to 
numerical methods. Here we will discuss the numerical method employed in this thesis. Let us con- 
sider a one- variable case; we can express the distribution function W(x, t) in terms of an appropriate 
set of basis function {p n (x)} 

W(x,t) = £V n (t)p n (s). (2.53) 

n 

The sum depends on the choice of basis function. Once we solve for the coefficients W n (t), we have 
full knowledge of the non-equilibrium distribution function. 



The use of Eq. (2.53) casts the Fokker-Planck equation into a set of recurrence relations 

Wn = -- - Qn'Wn-2 + QnWn-1 + Q n W n + Q+VF„ + 1 + Q+ + W n+2 • • • , (2.54) 

where Q's are some known constants, and we seek for short-ranged coupling. We will only encounter 
3-term recurrence relations in this thesis, namely 

Wn = QnW n -l + QnW n + Q+W n+1 . (2.55) 

This type of recurrence relations can be solved efficiently using the continued fraction method. The 
details of the continued fraction method are discussed in Appendix A.l. 
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Fokker-Planck versus Langevin 



Though it is possible to run Langevin simulations to obtain the average of a dynamical variable, 
solving Fokker-Planck equations with the continued fraction method requires much shorter com- 
putational time (few minutes on a laptop). The drawback is that we do not have any information 
about the trajectories. Though this drawback is offset by the fact that the distribution can also 
provide us valuable physical insights. 



2.5 Summary 

It is a long chapter, let us summarize what has been presented. We first described fluctuations 
and dissipation in an open system by modeling the bath as a set of harmonic oscillators. Using 
Hamiltonian mechanics, a Langevin-like equation of motion was obtained. Particularizing to the 
problems of particle and spin, we recovered the phenomenological Langevin equations. The example 
of Rubin Model provided us useful insights of the bath-of-oscillators model. 

We then derived the Fokker-Planck equation by making use of the continuity equation for the 
probability distribution. The Fokker-Planck equations for particle and spin were obtained from 
their corresponding Langevin equations. Eventually, the use of the continued fraction method in 
solving Fokker-Planck equations was discussed. In the next two chapters, we will demonstrate the 
use of this method in solving the Fokker-Planck equations for a particle in a periodic potential and 
a dipole (spin), both in the large damping limit. 
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Chapter 3 

Translational Brownian Motion: 
Particle in a Periodic Potential 

3.1 Introduction 

In this chapter, we apply the continued fraction method to solve the Fokker-Planck equation for 
a Brownian particle in a periodic potential. This problem finds applications in the non-linear 
pendulums, superionic conductors, phased-locked loops in radio, Josephson tunneling junctions, 
etc. |9]. Similar works can be found in [TU] and [TT], where the Langevin equation is used instead. 




_ 4 i 1 1 1 1 1 ] 

-15 -10 -5 5 10 15 

Figure 3.1: Periodic potential without (top) and with (bottom) biased force. 
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Consider a one-dimensional case, the particle is kicked around by the Langevin force. When the 
Langevin force is large enough, the particle will travel from one potential well to the next, causing it 
to diffuse in both directions. If we apply an external force, the particle will diffuse in one direction 
preferably (see Figure 3.2), and we are interested in the drift velocity (x). 




Figure 3.2: The trajectories of two independent Brownian particles in a periodic potential. The 
particles are subject to a constant force so that the random walk is biased to the force direction 
d 

The Langevin equation can be written as 

Mx + jMx + V(x) =F ext (t) + Z(t), (3.1) 
where F ext (t) is the external applied forc^J We consider a periodic potential of the form 

V(x) = -Focos(^x), (3.2) 
where 2Vq is the height of the well. The minus sign is inserted for convenience. 



3.2 Strong Damping: Smoluchowski Equation 

In the regime of high friction, the velocity of the particle reaches steady state rapidly, thus the 
inertial term, Mi, can be omitted. The resulting Langevin equation reads 

2-7T 2tt 

jMx(t) + — V sm(— x) = F exb (t) + f(t). (3.3) 

1 Risken [9] discusses the application to super-ionic conductors. A super-ionic conductor consists of a a nearly 
fixed ion lattice in which some other ions are highly mobile. If an external field is applied to a one-dimensional model, 



neglecting the ion-ion interaction, the equation of motion is the same as Eq. (3.1 1 
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Following Risken [9], we introduce the following dimensionless variables, 



2?r 



-x; t 



2vr V 



t; 7 



L /M 



7! ^c: 



L F Pr 



2tt Vr 



L L V M ' ' 2?r V Vb 

The Langevin equation is transformed to 

75 + sin(s) = F ext (t) + £(£), 

and the noise correlation function becomes 

(|(t)|(f)) = 27T<5(i-f). 







2?T Vn 



T 



Va ' 



(3.4) 



(3.5) 



(3.6) 



We will drop the tildes and it is understood that we are using the normalized units. The corre- 
sponding Fokker-Planck equation, according to (2.50), for the distribution W(x,t) is given by 



dW d 



7- 



dt dx 



sin(x) - F ext + T 



Ox 



W. 



(3.7) 



This is a special case of the Smoluchowski equation, the Fokker-Planck equation for an over-damped 
Brownian particle. 



3.3 Converting the Smoluchowski Equation into Recurrence Form 

In the steady state (long time limit), we expect the distribution to be periodic in space. In fact, in 
the problems of pendulums or Josephson junctions, the systems are indeed periodic (from to 2ir). 
Then, we can express the distribution function as a Fouries series in space 

00 

W(x,t)= J2 W n (t)e inx , (3.8) 



and we will need to solve for the coefficients W n (t). Substituting the expansion (3.8) into the 
Smoluchowski equation (3.7), we obtain the three-term recurrence relation 

7 W n = Q-Wn-l + QnWn + QiW n+ l, (3.9) 

where 

Q n = -inF ext (t) - n 2 T; (3.10) 

Qn = +2 n ; 

Qn = ~\ n - 

The zeroth term Wo, which is used as the "seed" in the continued fraction method (Appendix A.l), 
is fixed by the normalization condition. We normalize the distribution function over a period, 

W(x,t)dx = 1; W = ^-. (3.11) 



The drift velocity can be obtained by taking ensemble average of the Langevin equation (3.5), 

j(x) = F cxt (t)-(smx). (3.12) 
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Solving the recurrence relations, we get the distribution W(x, t) and hence the average of any func- 
tion involving x. From the expression above, we then can obtain the drift velocity (x). 

We will consider applied force of the form 

F cxt (t) = F + bcos(Qt). (3.13) 

The first term is a constant force which tilts the potential profile while the second term drives the 
system and allows us to study the dynamical properties. 
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3.4 Stationary Response (DC) 



In the absence of AC driving i ? ext (t) = F, the stationary solution (W n = 0) of the three-term 
recurrence relation (3.9 ) can be obtained using the scalar continued fraction method (Appendix A.l). 
In Figure 3^, the drift velocity is plotted against the applied force, F, at different temperatures. 
There are regions where the curves stay flat; the particle is "locked" in the potential well and there 
is not enough applied force or Langevin force (at low temperature) to push it away from the well. At 
large F (or high temperature), the effect of the potential well is less, and the velocity grows linearly 
with the applied force. In between, we have the depinning transition between the two regimes. At 
T ~ 0, this transition takes place when the force equals the cosine well depth (F = 1), and breaks 
the minima structure. 




Figure 3.3: Drift velocity against constant force F at various temperatures. The curves stay flat 
when the particle is trapped in the potential well. At large F or T, the curves are almost linear to 
F. 



3.5 System Under AC Driving 

In the presence of driving, the system will never reach a stationary state. Instead, in the long time 
dynamics, we expect it to be oscillating with the same period as the driving. This is similar to the 
case of a driven damped oscillator, where the driving frequency is the only time scale in the long 
time dynamics. Therefore, we can take care of the time dependence of the distribution function by 
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expanding it into a Fourier series in time {e lfcfM }. Together with the Fourier series in space {e znx }, 
we have 



oo oo 



W(x,t)= 2 W^ k) e inx e [km , 



n=~oo k=-oo 

In the driven system, we are interested in the susceptibility Xi defined as 

7<*>A = 7(*)(*)-7<*)o 



x (k) e +iknt + x *{k) e -iknt 



(3.14) 

(3.15) 
(3.16) 



where (x)o denotes the time-independent part of the drift velocity. We will adopt the following 
convention, 



x (k) = x >(k) _ ix "«. 



(3.17) 



In the regime of linear response (b — > 0), we only have the first order term (we omit the superscript), 
and the linear susceptibility is 



Xe +int + X*e~ int 



(3.18) 



The interpretation is easier in linear response: susceptibility is the coefficient of the time dependent 
part (of a dynamical quantity) that is oscillating at the driving frequency. The real part is in phase 
with the driving, while the imaginary part is the out of phase contribution. 



3.5.1 Perturbative Treatment 

In the case of weak driving, we can expand the distribution function in Taylor's series of the driving 
amplitude b, 



oo oo oo 



(3.19) 



q=0 n=-oo fc=-oo 

and we need to solve for the coefficients [q] . In this problem, the driving amplitude should be 
b < 0.2 for this expansion to hold. 



The expansion (3.19) turns the Smoluchowski equation ( |3.7| ) into 

QnWi%[q] + QM%] + Q+Wi%[q] = -fW[q], 

where 



Qn 
Qn 

Qn 

fn h) [Q] 



—inF — ikjQ — n T; 
1 

H — n\ 
2 



(3.20) 



(3.21) 



n 



w< L k+ %-i] + w< i k -%-i] 
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Figure 3.4: The linear susceptibility at T = 0.01. The curves of F = and F = 0.5 are barely 
distinguishable. The profiles change significantly for F > 0.9 as the particle is no longer trapped in 



the potential well and this corresponds to the sloped region of the curves in Figure 3.3 



The perturbative structure is clear in Eq. (3.20). The solution of the previous order equation [q — 1] 
enters into the inhomogeneous part of the next order [q]. The zeroth order equation is homogeneous, 
and its solution enters the right hand side of first order equation, and so on. This set of iterative 
equations can be computed up to any order, until the solution converges. 

We will first start looking at the linear response, i.e. we only solve up to the first order. The 
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linear susceptibility is plotted in Figure 3.4 Most of the interesting behaviors happen around the 
resonant frequency (the frequency of oscillations near the bottom of the wells, in our units, 7^ = 1). 
It is the same situation as the response of a driven damped oscillator. The increment of the constant 
field, up to F 0.9, does not change the susceptibility profile much; the particle is locked in the 
potential well for small fields. Note that we have used a low temperature T = 0.01, such that 



the locking behavior is evident (recall the deppining behaviors in Figure 3.3). When the field is 
increased until it reaches the depinning regime, the profile changes considerably. It is because the 
particle gains enough energy to hop from one well to the next, it is no longer trapped. 

Let us look at the low frequency region (adiabatic driving), where the system responds in phase 
with the field. At small force (up to F ~ 0.9), the response is zero. It can be explained by looking 



at the flat portion of the drift velocity curves in Figure 3.3 Turn on the driving, the system moves 
back and forth along the x-axis, but there is no vertical movement. At large force, the drift velocity 
curve is sloped, thus giving non-zero response. 



The real part of the susceptibility is proportional to the power dissipated into the medium due 
to the driving. The power is calculated by multiplication of the velocity and the force, P = F(t)v(t). 
In our context, 

b ^ 



P oc bcos(nt)^[ X e +[nt + x*e~ im 



(3.22) 



oc b 2 x cos 2 (fit) + x" sin(fii) eos(f^) 
Taking the time average over a cycle, one finds P oc x' '■ 



3.5.2 Exact Treatment 

Now we solve the Fokker-Planck equation exactly, using the matrix continued fraction method. This 
method is more computationally expensive, since it involves matrix multiplication and inversion. 



We define the vector 



wt 2) 
(+1) 



w, 



(+2) 



V 



/ 



Substituting the expansion (3.14) into the Smoluchowski equation (3.7), we obtain 



where 



W n _x + Q„W n + Q+ W n+ i = 0, 



inF - n 2 T)I - i^VLA - -nbB; 



Qn 


= (-inF 


Qn 




= 


Qn 





(3.23) 
(3.24) 
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and 



+1 



+2 



V 



/O 1 
1 1 

1 



B 



/ 







V 



1 

1 O 1 
1 / 



I is the identity matrix. The recurrence relations can be solved by the matrix continued fraction 
method (Appendix A.l). 

We plot the real time dynamical loops (hysteresis loops) : the curves of the drift velocity against 
the real time driving field, 6cos(Oi), in Figure 3.5 At small driving, the loops are elliptical since 
the response is linear to the driving field. When the driving field increases, the contribution from 



the (non-linear) higher harmonic susceptibilities, x 



(2) x (3) 5 x (4) 



etc., becomes significant and the 



loops are distorted. The area of the loop is proportional to the power dissipated. 

To see the contribution of the non-linear susceptibilities, we plot the first three harmonics in 
Figure 3.6 As opposed to the linear response, the susceptibilities at low frequency are lifted from 



zero at large driving field. It is because the constant and driving forces combined are large enough 
to kick the particle away from the potential well, the particle is no longer trapped. We also observe 
oscillatory behaviors of the susceptibility curves in the frequency range ~ 0.1 — 1. We will relate 
the oscillations with the celebrated Shapiro steps discussed below. 
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Figure 3.5: Top: Hysteresis loops at = 1 and T = 0.01 without biased force. At large driving 
force b, the elliptical shape is distorted. Bottom: The biased force Fq = 0.2 breaks the left-right 
symmetry. 
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First Harmonic, Real Part 



First Harmonic, Imaginary Part 




Second Harmonic, Real Part 




Third Harmonic, Real Part 
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Figure 3.6: The first three harmonics of the susceptibility, F = 0.5 and T = 0.01. Symbols are the 
results from the perturbative treatment. Oscillatory behaviors can be observed around ~ 0.1 — 1. 
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Figure 3.7: The time-average of the drift velocity at 7^ = 0.3 and T = 0.01. Shapiro steps occur 
at (x) = n^Q, and their widths increase with driving force. 



We plot the time averaged drift velocity (x) in Figure 3.7 The curves exhibit Shapiro steps at 
the multiples of the driving frequency, (x) = njQ, where n = 1,2 (see Ref. [T3] for a discus- 
sion). The system is locked at the resonant frequency = (x)/n. The appearance of the Shapiro 
steps is due to the fact that the quantity (x) becomes stable with respect to a small change in 
external parameter (F here) at resonant frequency [12] . One needs a finite change to push the sys- 
tem away from this stability. This phenomenon is well known in the context of Josephson junctions. 

To relate it with the oscillatory behavior observed in the susceptibility curves, we look at a 



particular value of F and vary the driving frequency, the curve in Figure 3.7 will rotate back and 
forth. The curve is flat when the frequency is resonant and sloped when it is not (see Figure 3.8). 
This "modulation" might explain the oscillatory features in the susceptibility curves. 



3.6 Summary 

In this chapter, we showed explicitly how the continued fraction method is used to solve a Fokker- 
Planck equation. We obtained both the time-independent and time-dependent solutions. We also 
presented two approaches of using the continued fraction methods: perturbative and exact. Strong 
non-linear effects are observed in the dynamical hysteresis loops under strong driving. We could 



include the inertia term M'x and solve the Klein- Kramers equation (2.51). This would require one 
more matrix index for the expansion in momentum, which entails larger computational efforts, see 
Refs. and [13]. 
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Figure 3.8: The time-average of the drift velocity versus frequency at b = 1, F = 0.5 and T = 0.01. 
One can relate the oscillations with the oscillations observed in the susceptibility curves in Figure 

EH 
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Chapter 4 

Rotational Brownian Motion: Debye 
Dipole 



4.1 Introduction 

Here we will study a Fokker-Planck equation of different structure, which involves a dipole. We 
investigate the problem of non-interacting dipoles subject to DC and/or AC fields. This problem 
was first studied by Peter Debye in the 1920's, and constitutes the first example of rotational 
Brownian motion. In the Debye model, he assumed high friction and isotropy, i.e. the Brownian 
motion exhibits no preferential direction. The orientation of the dipoles solely depends on the 
angle between the electric field and the dipole vector, 9 in Figure |4.1| This problem not only finds 
applications in dielectric relaxation, but also rotational relaxation of ferromagnetic nanoparticles 





P 



Figure 4.1: Pictorial skectch of a dipole, p in an electric field, E . 



4.2 Fokker-Planck Equation 



4.2.1 Debye Orientational Diffusion Equation 

In the high friction limit, the Fokker-Planck equation describing the dipoles is [2] 




(4.1) 
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where £ is the viscosity coefficient. We define the Debye relaxation time and the dimensionless field 
parameter 



TD 



c 



a(t) 



pE(t) 



2k B T' ~~ y ~' k B T 
Making the transformation x = cosO, the resulting Fokker-Planck equation reads 



2t d 



dW(x,t) 



d_ 
Ox 



1-x' 



o 



a(t)\ W(x,t). 



(4.2) 



(4.3) 



4.2.2 Method of Solution 



Because of the range of x (—1 to +1), the natural choice of the basis function here is the Legendre 
polynomials P n (x)R 



W(x,t) = J2 W n(t)Pn(x). 



n=0 



Substituting Eq. (4.4) into Eq. (4.3), we obtain the three-term recurrence relation 



2 TD W n = Q~W n -l + Q n W n + Q+ +1 W n+1 , 



where 



Qn 
Qn 

Qn 



n(n + 1); 
n(n + 1) 



+a{t) 
-a(t) 



2n- 1 ' 
n(n + 1) 



(4.4) 

(4.5) 

(4.6) 
(4.7) 



2n + 3 

Similar to the previous chapter, we consider the combination of DC and AC fields, 

E(t) = Eq + E d cos(nt); 

a(t) = a + bcos{iU); a = T — ; b=——. 

k B T k B T 

The object of interest here is the average orientation (x) = (cos 6). 



(4.8) 
(4.9) 



4.3 Equilibrium Properties 



Without driving, we solve for the equilibrium solution to the recurrence relation (4.5) with the 



scalar continued fraction method (Appendix A.l), and study its average orientation. In fact, we 
1 The Legendre polynomials can be expressed as the Rodrigues' formula 

They obey the orthogonality relation 

f+l P n (x)P m (x)dx = 2^pr<5" 
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can obtain the analytic result using elementary statistical mechanics. The canonical distribution of 
such a system is given by 



-0H 



z 



(4.10) 



where the Hamiltonian is H = — p • E and Z the partition function^} The average orientation is 
defined as 

(cos 9) = / dT Pc cos9, (4.11) 



(4.12) 



where the integration is carried over the solid angle. The result is the Langevin function 

1 

ao 



(cos 9) = L(ao) = coth(ao) 



The Langevin function is plotted together with the numerical results in Figure 4.2 At small 

field, the average orientation grows linearly with ao since L(«o) ~ a o/3 — afj/45 + At large 

field, we reach the saturated region where (cos 9) ~ 1, the dipoles are almost fully aligned with the 
field. The comparison with the Langevin function serves as a test for our numerical method. 




Figure 4.2: The average orientation of a dipole, the solid line is the numerical results while the 



symbols are from the Langevin function Eq. (4.12). 



One can easily show that Eq. (4.101 is the stationary solution by substituting it into the Fokker-Planck equation 



(4.11. 
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4.4 Driven Dipole 



As in the previous chapter, we focus on the time periodic solution and perform the Fourier time 
expansion of the distribution function of a driven system, 



oo oo 



W{x,t) = Y^ Yl W^e iknt P n (x). (4.13) 

n=0 k= — oo 

We will not discuss the perturbative treatment and the exact treatment again, they are almost 
identical as in Chapter 3. 



Again, the susceptibilities \ are defined as 

(x)a = (x(t)} - (x) 



£ 

k=l 



x (k) e +iknt + *(k) e -ikcu 



(4.14) 



where (x)q is the time-independent part. In the regime of linear response (weak driving), we only 
keep the linear susceptibility 

b 



(x) A 



Xe +[nt + x* e~ ikQt 



(4.15) 



4.4.1 Linear Response 



From linear response theory, the susceptibility at zero DC field is given by the Debye relaxation 
formula 



1 



1 



X 



3 l + iftr D 



(4.16) 



The expression is compared with the numerical results in Figure 4.3 as a check of our numerical 
method, the top curves (ao = 0) of both panels. We also try a heuristic expression for the linear 
susceptibility at non-zero DC field, 

X = L'(ao) - ■ * , (4.17) 



1 + i f2 T cS 



where the effective relaxation time is 



TeS = T D 



L(a ) 



1 



«0 



-L(a ) - L 2 (a ) 



(4.18) 



(The effective time is defined as the initial slope of the relaxation curve when there is a small change 
of the applied field [TB].) 

We compare the expression with the numerical results in Figure |4.3| and good agreement is 
observed. The reason the derivative of the Langevin function, L'(ao), enters can be justified. 
Without driving, the average orientation takes the value of L{olq). Turn on the driving, we move 
back and forth along the Langevin curve and the vertical movement depends on the slope L'(cto) 
of the curve. = 0) decreases with increasing constant field because of the decreasing slope of 

the Langevin function. In the limit when the dipoles are nearly fully aligned by the large constant 
field, it is more difficult to rotate them, and their response drops. 

Note that the power dissipated is proportional to the imaginary part of the susceptibility, as 
opposed to the real part in the transport problem. It is due to the fact that we need to take the 
time derivative of {x)a to get the "velocity", giving (x)a oc if2 {x)a- The imaginary unit exchanges 
the role of the real and imaginary parts in the dissipation. 
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4.4.2 Beyond Linear Response 



Beyond linear response, we solve for the polarization at arbitrary AC field, and the hysteresis loops 



are plotted in Figure 4.4 The deformation of the elliptical loops can clearly be seen at large driving, 
due to the contribution of higher harmonic susceptibilities. The loops develop a spike- like structure 
at the tips, as in the custom hysteresis loops of magnetism. 

We also plot the first three harmonics of the susceptibility in Figure |4.5[ The susceptibility can 
no longer be described by the phenomenological Eq. (4.17) at large driving (it only works for b — )■ 0), 



as one can see from the susceptibility curves. As we increase the driving amplitude, the magnitude 
of the susceptibility at low frequency (adiabatic) drops. This is also related to the decreasing slope 
of the Langevin function. 



The results from the perturbative treatment are plotted in Figure 4.5 as a consistency check. 



Real part 




Figure 4.3: The linear susceptibility. The symbols are from the expression (4.17) while the solid 



lines are the numerical results. Both real and imaginary parts decrease at large constant field, ao- 
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0.8 




b cos(Q t) 

Figure 4.4: Hysteresis loops at t^Q = 0.5 without (top) and with (bottom) constant field. The 
elliptical loops are deformed at large driving field, b = f^- In the bottom figure, the constant field, 
oq = 0.5, breaks the left-right symmetry. 



The radius of convergence for the perturbative treatment is about b ~ 3, which is large as compared 
to the previous chapter where b ~ 0.1. In this range, we can already observe significant deviation 
from the linear response results (b = 0) before it breaks down (symbols in Figure 4.5). In the 
problem of Brownian particle in a periodic potential, the perturbative expansion fails before we can 
observe any significant non-linear effect. 
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Third Harmonic, Real Part 




Third Harmonic, Imaginary Part 




Figure 4.5: The first three harmonics of the susceptibility at oto = 1. The symbols are the results 
from the perturbative treatment, which fails at b > 3. The curves of b = represent the linear 
response results, and are related to the derivative of the Langevin function. 
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4.5 Discussion 



Here we end our discussion on classical open systems. To conclude, we have presented the Hamil- 
tonian viewpoint of open systems and two equivalent ways of studying them: the Langevin and the 
Fokker-Planck equations. We then used the continued fraction method to solve the Fokker-Planck 
equations for particle in a periodic potential and dipole under driving force, and studied their linear 
and non-linear responses. We shall devote the next two chapters to the study of open quantum 
systems. 



35 



Chapter 5 

Quantum Open Systems 



5.1 Hamiltonian and Reduced Description 

In quantum open systems, we study the fluctuation and dissipation on a quantum system due to the 
interaction with the environment. Unlike the classical counterparts, it is difficult to introduce phe- 
nomenological equations to describe these effects, because of the unitarity of the quantum dynamics. 
Previous attempts are plagued with various problems, e.g. violating the uncertainty principle or 
the superposition principle [5]. 

Therefore, it is natural to view an open quantum system as a system with a few degrees of 
freedom coupled to a bath with many (infinite) degrees of freedom. The total Hamiltonian, H to t, 
describing such a model is written as 

Hot = H sys ®I + I®H h! , th + V iQt (5.1) 
= H + V int , 

where H is the free Hamiltonian of the system plus bath and Vint is their interaction Hamiltonian. 
The combined system is fully described by the total density matrix, ptot > whose dynamics is governed 
by the von Neumann equation 

^ = (5.2) 



However, in most cases, we are only interested in the properties and evolution of the system 
(we do not have much control on the bath). We then introduce the reduced density matrix of the 
system, obtained by partial tracing over the bath degrees of freedom, 

p = Tr B [ptot]- (5.3) 

The reduced density matrix contains all the information we need in most cases of interest (heat 
transport is an exception though). In the following sections, we will derive a differential equation 
to describe its time evolution using the second-order perturbation theory. 
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5.2 Perturbation Theory in System-Bath Coupling 
5.2.1 Evolution Operator 

In general, we are not able to obtain an exact equation of motion for the reduced density matrix 
similar to the classical equation Eq. (2.17); we have to resort to perturbation theory. We start with 
the evolution operator, 



e — %Htat (*— *o)_ 



U(t,t ) 

To facilitate the perturbative treatment, we use the identity [16 



e P{A+B) _ Q I3A 



1+ / d Xe~ XA Be^ A+B ^ 



(5.4) 



(5.5) 



It can be confirmed by multiplying both sides by e ^ and differentiating with respect to (3. Using 
this identity, the evolution operator can be expressed in the following Dyson-like form 

rt-to 



U(t,t Q ) 



-jiH(t-t ) 



1 



dse L n Hs V^ t e-^ H+v ^ s 



(5.6) 



Let us make the transformation s — > s — to and keep up to the first order term in Vint) we obtain 



£/(Mo) « U (t,t ) 



ds Vint(s) 



to 



where 



U (t,t ) 



rt 



U J Q (s,to)Vi n tU (s,to). 



(5.7) 



(5.8) 
(5.9) 



5.2.2 Heisenberg Equation 



Now instead of looking at the evolution of the reduced density matrix, we look at the Heisenberg 
equation of an operator acting only on the system's Hilbert space, 



dA(t) 
dt 



-[H tot ,A(t)} 

^[F sys (t),A(t)] + ^[M nt (t),A(t)]. 



(5.10) 



Other than telling us the evolution of an operator (say momentum or position operator of the sys- 
tem), the above equation also allows us to study the evolution of the reduced density matrix upon 
choosing an appropriate operator, as what we will do in the next section. 



Making use of Eq. (5.7) and the similarity transformation property 

[V int (t),A(t)} = 0% t ) [V^t, A]U{t,t ), 



(5.11) 



the Heisenberg equation becomes 

dA(t) 
dt 



i[H S2/s (i),^(i)] + ^V in t(i),I(i)] 



j^ ds { Vnt (S) [Vint (*M(t)] + [A(t), Vint (*)] Vint (*) }• (5-12) 



37 



The second order perturbative structure is clear. The first term arises from the free evolution, the 
second and third terms are the first and second order corrections, respectively. It is worth recalling 
that we are still in the Heisenberg picture, while the operators with tildes are evolved by its free 



Hamiltonian [cf. Eq. (5.9)] 



We have obtained a generic equation of motion for a system operator based on perturbation 
theory. In the next few sections, we will make use of this equation to derive the so-called master 
equation: the equation of motion for the reduced density matrix. 

5.3 Quantum Master Equation (Bloch-Redfield) 
5.3.1 Hubbard Operators 

We introduce the Hubbard operators X nm = \n)(m\, where {\n}} is the set of energy-eigenstates of 
the system Hamiltonian, 

H sys \n) = e n \n). (5.13) 

The properties of the Hubbard operators can be found in Appendix A. 2. The Heisenberg equation 
of the Hubbard operator is 

where the energy level difference is A nm = e n — e m and the relaxation term reads 

Rum = -pjf ds{^ nt (s)[^ nt (t),X nm (t)] + [X nm (t),^ nt (t)]^ nt (s)}. (5.14) 

We introduced the Hubbard operators because they are useful in getting the elements of the 
reduced density matrix by tracing, 

Pmn{t) = Tr s p(t )X nm (t) . (5.15) 

Before we could obtain the master equation, we need to make a few more assumptions. 

Decoupled Initial Condition 

The factorized initial condition 

Ptot(*o) = P(t ) <8> Pbath(*o), (5-16) 

is a vital assumption in the process of tracing and getting the matrix elements. We will discuss the 
decoupled initial condition at the end of this chapter. 

Coupling Structure and Bath "Centering" 

The coupling is in standard factorized form (additional summation would give the general case, but 
this just brings notational changes), 

V mt = F®B. (5.17) 

The coupling is defined such that the first moment of the bath Hamiltonian vanishes, (B) = 0. If 
this were not the case, we redefine the coupling, and lump the non-zero average into the system 



Hamiltonian [this is equivalent to the zero average of the Langevin force in Eq. (2.21)]. 
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5.3.2 Bloch-Redfield Equation 

With these conditions, the equation of motion obtained after tracing is 



dt 



where the relaxation term is 



Rr, 



li 2 



^ Pm'n'(t) 



mm) F n ' n F mm i 



(W* n , + w, 

- ^ 5mm'Fn'lFl n Wi n , 

I 

- ^2 ^nn'FmlFlm'Wim 



(5.18) 



(5.19) 



in which (n\F\ 771) — F nrn . Th.6 transition rate is 

W nm = f d se iA„m( S -t)A (B(t)B(s)) 

Jtn 



f ^ dTe-' lA ™ T/h K(T), 
Jo 



(5.20) 
(5.21) 



where the bath correlator is K[t) = (B(t)B). The real part of the transition rate determines the 
speed at which the stationary solution is reached. Upon taking the initial time at to — > — oo, the 
transition rate becomes a half Fourier transfrom of the bath correlator. 



Equation (5.18) is frequently called the Bloch-Redfield equation [17J, which is widely used in 



magnetic resonance (nuclear, electron, or ferromagnetic), optical spectroscopy, laser physics, and 
electron-transfer reactions in molecules and bio-molecules. This equation is generic, we have not 
specified our system, bath or the coupling. All the information of the bath goes into the correlator 
K(t), and we will need to specify the bath coupling, B, the spectral density J(uj) and its initial 
state Pbath^o)- 

5.3.3 Ladder Couplings 

Let us consider the couplings of the type 

(n\F\m) = F nm = L~5„ im _i + L+(5 n/m+ i. (5.22) 
It is suited to study the harmonic oscillator problem with coordinate/momentum coupling 

F = r ? (a + a t ), (5.23) 
or spin problems with general couplings of the form 

F = v+ {v(S g ),S-}+ri-{v(S g ),S + }, (5.24) 
where {•, •} is the anti-commutator and rj± describes the symmetry of the interaction [21J. 
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With such coupling, the relaxation term becomes 
fi 2 Rmn = — L n+1 W* +l n + L^ l _ 1 L n W*_ ln + L^L m+1 W m+ i jm + I/^__ 1 L m H / m _i jTra ^ p mn 

+ (Wn,n-l + Wm,m-l)Ln L m-\ Pm-l,n-l 

,71+1 + ^m,m+lj L m+1 p m +l,n+l 

,71-1 ,m+l) ^m+1 Pm+l.n-l 

+ (Wra,ra+l + ^m,m-l) -^m-1 Pm-l,n+l 

~W*_ l n _ 2 L n _ 1 L n Pm,n-2 ~ ^r!+l,7i+2^n ^n+1 Pm,n+2 

-W m -i. m - 2 L^ n _ 1 L^ n _ 2 p m -2,n ~ W m+ i. m+2 L m+1 L m+2 p m +2,n- (5.25) 

5.3.4 Secular Approximation 

At this stage, one invariably invokes the secular approximation, discarding the terms of the type 
L + L + and L~ L~ . We basically get rid of the last four lines of the relaxation term above and keep 

fr 2 Rmn = ~(L+L- +1 W* +ltn + L^_ 1 L~W*_ l n + L+L~ l+1 W m+1 , m + L+_ 1 L~ W m _i, m ) p mn 

,771-1^-^71 ^777-1 Pm— 1,77— 1 
,71+1 + ^m,m+l) "^77 -^777+1 P771 + l,rt+l- (5.26) 

We will justify this approximation in the discussions at the end of the chapter. After the secular 
approximation, the coupling of the matrix elements becomes simpler. Any matrix element, p mn , is 
only coupled to its adjacent diagonal neighbors, p m+ i in+ i and p m -i,n-i- This short-ranged coupling 
simplifies the implementation of the continued fraction method. 

5.4 Application to the Bath-of-Oscillators Model 
5.4.1 Hamiltonian Redux 

In this model, the system is coupled linearly to the coordinates of a bath of oscillators (B = 
J2a c a x a), as in the classical case. The Hamiltonian is 



Htot — H sys 



N 2 2 

^— ' V2m a 2 V m a w; / J 



Because of the close resemblance to the classical counterpart, this model gains the name of "Quan- 
tum Brownian Motion". The only difference is that all x's and p's are now quantum operators. For 
instance, we will substitute the bath coordinates with the bosonic operators, x a oc ^ ^aj, + . 

5.4.2 Bath Correlator 

The bath is assumed to be at thermal equilibrium initially, 

g— P -Hbath 

Pbath(to) = 79 • (5-28) 
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The bath correlator thus reads 
K(t) = (B(t) B) 



where the Bose function is 



h E W^^^J^ + a a e^ T )(ai + a a )) 

a " 
poo J 

ft / — J(w) n w e iwr + K + l)e-^ 

f°° doj r 1 
h l — Ji^) cothf— /3fto;) cos(o;r) — i sin(a;r) 
Jo 7T L 2 



= l/(e^ - 1). 



(5.29) 



(5.30) 



At high temperature, the real part of the correlator behaves classically, as in Eq. (2.25). The 



imaginary part arises from the non-commutability of the bosonic operators and has no classical 
analog. 



5.4.3 Relaxation Coefficients 
Counter- Term and Renormalized Rate 

The counter-term F 2 , can be handled by modifying the transition rate as 

W ^ W + ih^(0), 



where 7(0) is the damping kernel evaluated at t = [cf. Eq. (2.26)] 

duo J(co) 



7(0) 



1 

M. /n 



7T CO 



(5.31) 



(5.32) 



Drude-Ohmic Damping 



We discussed in the Rubin model (cf. Section 2.3.4) that, in a physical model, the density of the 
bath modes is cut-off at high frequency. A frequently used model satisfying this condition is the 
Ohmic spectral density with Drude cut-off, 



J{u) 



1+oj 2 /uj 2 d 



2 ' 



(5.33) 



The damping coefficient, 7, keeps track of the coupling strength, and is proportional to c 2 in the 



total Hamiltonian (5.27). This should not be confused with the gamma with time argument, which 



is the damping kernel. 

For the above spectral density, transition rate is found to be 
W mn 
Re[W(A)] 

J. t iujjjj )~ c — — l 

7 A 2 7A f ,,PfruD 



Im[W(A)] 



W(A mn ); 

7A 1 
1 + (A/huj D ) 2 e^-l 
A 2 



(5.34) 



Rc 

2huo D l + (A/hw D ) 2 Tr[l + (A/hw D ) 2 ] l rv 2tt ' Wy 2v' n ' f3hw D 



+ 



i\) is the Psi function defined as the derivative of the logarithm of gamma function, if)(x) = Jdnr(x). 
Some properties of the transition rate are discussed in Appendix A. 3. 
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5.5 Discussion of the Approximations 



At this point, we have our master equation ready for use, we just need to specify the spectrum 
(energy levels) of the system. But before we try to solve for any system, there is a need to discuss 
some of the subtleties of the master equation. 



Weak Coupling 

We need to be more specific on what we mean by weak coupling. We assume Vint is small, and thus 
R mn can be treated as a perturbation to the free evolution. But letting \t — to\ — > oo, the integral in 
the transition rate W [see Eq. (5.21)] would become very large and the perturbation theory breaks 
down. In many problems of interest, there exists a correlation time r c , such that the correlator is 
negligible K{t c ) ~ after r c . Thus, the integral would not grow as we feared. The weak coupling 
approximation is then valid in the regime of 

7T C «1. (5.35) 



Secular Approximation 

The secular approximation is equivalent to the rotating wave approximation (RWA) in quantum 
optics. Specifically, a coupling of the form, V oc F^(a+a))+F + (a+a)) is reduced to V oc F_a) +F + a. 
It is argued that the secular term F + a) is rotating at Q 1 ( UJ +^) t ^ which is much faster than the term 
F^a) oc e 1 ^^)*, and can be averaged out. In the case of A being negative, one discards the terms 
F-(J and F + a instead, since they are the ones which are rotating faster. 



Decoupled Initial Condition 

We have chosen the factorized initial condition to facilitate the process of tracing. In the context of 
condensed matter, the system and bath do not meet at our chosen time, they have always been in 
contact. Thus, we set the initial time at minus infinity, and hope that the artificial initial condition 
is forgotten by the time we start to manipulate the system at t = [5]. 

In the next chapter, we will use the master equation to study a damped quantum harmonic 
oscillator and make comparison with some exact results as bench-marking. 
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Chapter 6 

Quantum Harmonic Oscillator 



6.1 Introduction 

In this chapter, we study the properties of a quantum harmonic oscillator linearly coupled in coor- 
dinate to the bath. The Hamiltonian is 

2 1 N 2 1 2 

"... = £ + i"*? + E[|t + " sSg*) ]■ 

Q=l 

The study of the damped quantum harmonic oscillator is important because this model applies to 
any system slightly displaced from its stable local potential minimum. On the other hand, this is 
one of the few problems amenable to exact solution, by the use of path integrals or diagonalization 
of the Hamiltonian. Thus, we can make comparison and assess the validity of the master equation 
derived in the previous chapter. This comparison is essential before we solve for systems with no 
exact solution. 

We will first show how to cast the master equation into a set of recurrence relations and how to 
solve them. After which we will study the equilibrium and driven properties of a damped quantum 
harmonic oscillator. 

6.2 Method of Solution 

6.2.1 Coefficients of the Master Equation 

The eigen-energies of a quantum harmonic oscillator are equally spaced, 

e n = (n + 77)^0; A nm = (n - m)hujo. (6.2) 



The matrix elements of the coupling Hamiltonian (F = x) entering the relaxation term (5.26) are 



n 
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6.2.2 Casting the Master Equation into Recurrence Form 

We can construct vectors from the columns of the reduced density matrix (truncated at the N th 
level) , 

/ POn \ 

Pin 
P2n 



\ PN,n J 

and the master equation acquires the following recurrence form 

Qn c n-l Qn c n + Qn c n+l 

where the elements of the matrices Q's are 



(6.5) 



(Qn)m,m 
(Qn )m,m+l 



n 



-A 



1 



-(w* ,1 

-2 y rr n,n+l 



(Qn 



m,m—l 



ft 2 
1 



+ W m>m +i )I/ n L m+1 ; 



W n,n-l + W m ,m-l )L n L m _ v 



6.2.3 Implementation 

The truncation level N has to be chosen such that NhuJo/ksT 3> 1, and the energy levels higher 
than N become irrelevant. When we solve for systems with finite levels (e.g. spin problems |21j). 
there is no need to perform truncation and the recurrence relation is exact. 

We are happy when we see the 3-term reccurrence relation, since we know how to solve it with 
the continued fraction method. To obtain the solution, we first need the "seed" Co so that all other 
c n 's can be calculated by the equation (cf. Appendix A.l) 



n L n-l 



+ a r , 



(6.7) 



But there is a problem in solving for the stationary solution (c n = 0). Eq. (6.5) is a set of 
homogeneous equations where the solution involves a multiplicative constant. Unlike the classical 
problems, we cannot obtain the "seed" Co, from the normalization condition Tr[p] = 1, as it involves 
all the vectors c n 's. This problem can circumvented by fixing one of the matrix elements in Co- The 
first vector obeys 



Ac = (Q + Q Si)c = 0, 



(6.8) 



where Si is some function of the Q's. We provide an extra equation by requiring the first element 
of Co to be 1, getting a set of over-determined equations 

A 
d 



co 



where 



(i, o, o, o; 









V i J 
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We can solve this set of equations by using the least square method. The subsequent c n 's can then 
be generated by the upward iterations Eq. (6.7) in the continued fraction method. Eventually, we 
will need to normalize the density matrix by the operation p = p/Ti[p\. 



As opposed to the classical problems, the indexed structure is automatically given by the energy 
levels, it saves us the troubles of choosing an appropriate basis function and manipulating to get the 
recurrence structure. The price to pay is the extra effort in getting the seed. The master equation 
can actually be solved by inverting a matrix of dimension A'" 2 x iV 2 , which involves computational 
efforts of 0(iV 6 ). In using the continued fraction method, all the matrices are of dimension N x N, 
thus we have reduced the complexity to 0(A 3 ) with iV iterations. This allows us to reach the high 
temperature regime, where high energy levels are excited and iV becomes large. In the next few 
sections, we will use the continued fraction method to study the equilibrium properties and the 
response to applied fields. 



6.3 Equilibrium Results: Dispersion of Coordinate and Momen- 
tum 

The equilibrium solution (c n = 0) is found to be the canonical distribution 

g— /3 -ffsys 

Peq = = , (6.9) 



independent of the coupling strength Some of the interesting thermal-equilibrium quantities 
are the mean square of the coordinate and momentum, plotted in Figure |6.1| The results from the 
canonical distribution are 

^V> = jj-(p 2 ) = coth(/3^ /2). (6.10) 

At high temperature, both quantities grow linearly with the temperature. The system behaves 
classically at high temperature and follows the equipartition law. At T = 0, they reach the standard 
quantum limit \J (x 2 )(p 2 ) = h/2, as opposed to zero in the classical case, a manifestation of the 
zero point energy. 

Therefore, we have achieved our goal of reproducing the dispersion curves by solving the quantum 
master equation with the continued fraction method. The results are in full agreement with the 
statistical mechanics, using a truncation level oi N = 30. 



6.4 Response to Static Field 

The master equation after the secular approximation is of Pauli type: it involves only the diagonal 
terms. In order to check our handling of the off-diagonal structure, we apply a static force, Fq, to 
the harmonic oscillator. We assume the applied force is small (of the order of 7), such that we can 
ignore the change to the relaxation term. R mn is already of the order of 7, thus, any modification 

1 Though the numerical solution does not depend on the coupling strength, the exact solution says otherwise [5]: it 
is canonical only in the limit of 7 — > 0. For any finite 7, the exact solution is different from the canonical distribution, 
and the correction is of the order of 7 |18| . To explain the discrepancy, one should recall the approximations we have 
made: the secular approximation and the weak coupling approximation. Because of these two approximations, we 
always obtain the canonical distribution as the stationary solution. However, the correction is small within the weak 
coupling regime. 
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Figure 6.1: The dispersion curves of the coordinate and momentum. The numerical results coincide 
with the canonical distribution. At high temperature, both quantities approach the classical limit. 



is at least of the order of j 2 and can be discarded. We just have to alter the free Hamiltonian, the 
equation of motion of the Hubbard operator then reads 



dXr, 



dt h 

The corresponding master equation is 

dpmn i 



H-xF ,X n 



+ Rr. 



dt 



^nmPmn ~i~ D mn + R m n > 



where R mn remains the same as Eq. (5.26|) and 

iF 



D r , 



(m+ 1) 



Pm+l,n 



[n + 1) 



~ Pm,n+1 



Pm-1, 



Pm,n—1 



(6.11) 



(6.12) 



(6.13) 
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gives the off-diagonal structure. The addition of the DC force does not break the short-ranged cou- 
pling of the matrix elements, we are still able to solve Eq. (6.12 ) with the continued fraction method. 



The DC response is characterized by the DC susceptibility defined as 

Xdc = (x)/F . (6.14) 

Let us first see how this quantity behaves classically. The Langevin equation of a forced oscillator 
is 

x + 7± = £(t) +F - MujI x. 
At equilibrium, both terms on the left hand side vanish. Taking average, one finds 



1 



Xdc 



(6.15) 



(6.16) 



which is independent of the temperature. As the system is linear, we expect this to hold in the 



quantum problem. Figure 6.2 shows that this is indeed the case, and it serves as a check of our 
handling of the off-diagonal structure of the quantum master equation. 
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Figure 6.2: The DC susceptibility at 
independent of the temperature. 



F 



0.0001, 7/w = 0.01 and oj d /u = 10. It is 



6.5 Linear Response to Time-Dependent Field 
6.5.1 Perturbative Chain of Equations 

As we did in the classical problems, we apply a small AC field, FdCos(Qt), and study the system's 
response characterized by the AC susceptibility 

(x) = Y{xa C e +int + x: c e- int ). (6.17) 
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To obtain the new master equation, we just need to replace Fq in Eq. (6.11) and Eq. (6.12) by 
FdCOs(p,t). Assuming the driving force is small, the system is only slightly perturbed from its 
equilibrium state. We can split the matrix elements into time-independent and time-dependent 
parts 

o - JO) + E± ( (+ 1 ) e +int + o (_1) e- in ^ (6 18) 

Pmn — Pmn ' \"mn e 'Pmn e Ji \v.LO) 

and solve the recurrence relations for zeroth order and first order in sequentially. 

Zeroth Order: ^A nm + R^ n = 0, (6.19) 

First Order: if A nm /h - + R^ n = -D^jF d . (6.20) 



The superscript in R mn and D mn indicates the order of density matrix element that should be used 
in the expression. The solution of the zeroth order enters into the first order equation, on its right 
hand side. Therefore, these equations have to be solved sequentially. 

6.5.2 AC Susceptibility Curves (Dispersion and Absorption) 

The AC susceptibility is plotted in Figure [673] and we observe resonant behavior at £1 ~ ujq. The 
imaginary (absorptive) part is of Lorentzian shape with peak at around ujq. At resonant frequency, 
the real part becomes zero because the system is oscillating completely out of phase with the driving 



((x) oc sin fit), as one can see from Eq. (6.17) by setting x' = 0. The dependence of the absorptive 
part on the damping strength can be seen in Figure [674| The peaks become flattened and shifted to 
the right when the damping strength is increased. This dependence is useful if we wish to estimate 
the damping strength by measuring the height or the full width half maximum of the absorption 
peaks. 
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Figure 6.3: The AC susceptibility at = 1, 7/wo = 0.01 and ujd/ujq = 10. Resonant behavior 
occurs around Q/coq 1. 
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These results are in good agreement with the analytical result obtained by stochastic modeling 

® 

11 7 
X=^- T — 2 2 : — ^7"^' 7M = i = — 1 ■ ( 6 - 21 ) 

Since the agreement is nearly perfect, we have not shown the analytical curves in the figures above. 



y/™ =0.01 
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Figure 6.4: The imaginary part of the AC susceptibility at = 1 and ujd/ujq = 10. The 
absorption peaks flatten when the damping is increased. 



6.5.3 Secular Approximation Revisited 

As one can see from the susceptibility curves, the dominant contribution comes from the region 
around the resonant frequency. We can relate it with the secular approximation we have made 
earlier. We keep the terms with resonant frequency and discard the off-resonant terms since their 
contribution is minute. As one increases the damping, the curves flatten and the contribution from 
the region away from the resonant frequency becomes important as well. This means the secular 
approximation is less accurate at large damping. Actually, the error introduced by the secular 
approximation is of the order of 7 2 [19]. A thorough investigation of the secular approximation in 
the Jaynes-Cummings model (two-level system in one-mode electromagnetic field) in a dissipative 
bath can be found in Ref. |20j . 

6.6 Discussion 

To conclude, we have used the master equation of Chapter 5 to solve for the damped quantum 
harmonic oscillator and discussed some of the implementation issues. We studied its thermody- 
namical properties with and without constant field. After which we investigated the response to 
a time-dependent field and obtained the absorption-dispersion curves. The results are consistent 
with the exact results, provided 7/o;q <C 1. 
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This agreement gives us confidence on our handling of the quantum master equation and our 
implementation of the continued fraction method. It is important as this efficient method is inde- 
pendent of the approximations used to get the quantum master equation, so our implementation 
could be extended to master equation with improved approximations. Possible extensions include 
the removal of the secular approximation or the inclusion of higher order terms in the system-bath 
coupling. These will break the 3-term recurrence structure, and more efforts are needed to solve 
them. The investigation is still ongoing and we shall leave them out of the thesis. 
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Chapter 7 



Summary 



To end the thesis, let us summarize our main results: 

• In the classical regime, we reviewed how the fluctuation and dissipation of an open system can 
be explained by modeling the bath as a set of harmonic oscillators. Two equivalent approaches 
were presented to study classical open systems: the trajectory approach (Langevin equation) 
and the distribution function approach (Fokker-Planck equation). The use of the continued 
fraction method in solving Fokker-Planck equations was also discussed. 

• We then demonstrated the use of the continued fraction method in solving the Fokker-Planck 
equations for particle in a periodic potential and classical dipole, both in high friction limit. 
We solved for both equilibrium and time-dependent solutions. The time dependent solution is 
obtained at arbitrary DC and AC fields. At large AC field, we observed significant deviation 
from the linear response results. The non-linear effects are reflected in the distortion of the 
shapes of the dynamical hysteresis loops. 

Previous studies of the particle problem focused on the zeroth and first harmonic suscepti- 
bilities; and for the dipole the first few harmonics but in the limit of weak driving. Here we 
obtained all the harmonics, for any driving strength and biasing DC field. We interpreted their 
characteristic features and presented all the information in a compact way with the dynam- 
ical hysteresis loops; this was not discussed in the literatures before. The assessment of the 
perturbative approach versus exact approach is also new for both problems. This assessment 
can be of valuable methodological interest. 

• In the quantum regime, we reviewed the derivation of the quantum master equation for the 
reduced density matrix, and its application to the bath-of-oscillators model. Approximations 
and subtleties of the master equations were also discussed. 

• We went on to solve the master equation of a damped quantum harmonic oscillator using the 
continued fraction method. This problem allows us to make comparison with exact results 
and assess the validity of the master equation. We investigated both the equilibrium and 
time-dependent solutions. Driven systems are more challenging in the quantum case, so we 
are content with implementing a linear response treatment. These results are in good agree- 
ment with the exact results (when available) under the condition 7/^0 *C 1. 

This showed that we successfully implemented the efficient continued fraction method in 
solving the master equation. This method reduces the computational complexity significantly 
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and allows us to solve for systems with many levels. In fact, the problem of damped quantum 
harmonic oscillator constitutes the first attempt in using the continued fraction method to 
solve for a mechanical system (it was originally proposed and tested for spin systems with 
finite number of levels). This widens the application range of this approach, and opens doors 
for further studies. 
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Appendix A 

Appendices 



A.l Continued-Fraction Method 

Here we give a summary on solving the 3-term recurrence relation of the form 

QnCn-i + Q n c n + Q+c n+ i = -/„ , n = 0,1,2,3, , (A.l) 

with the continued fraction method. The coefficients Q's and the inhomogeneous part /'s are some 
known constants. Risken [9] introduced the following ansatz 

c n = S n c n —\ + a n , (A-2) 

and obtained the relations 

g On . a _ Qn a n+l + fn g\ 

Qn "f" Qn £>n+l Qn ~\~ Qn $n+l 

For finite recurrence, c n >;v = for some N. We can enforce this by setting = 0, ajy = 0, and 



generate all the other 5 n 's and o n 's by the downward iteration Eq. (A. 3). To obtain all c n 's, we 
only need the "seed" Cq, which can be obtained by solving 

(Qo + QtS 1 )c = -(fo + Q^a 1 ). (A.4) 



Other c n 's can then be generated by the relation Eq. (A. 2). In the case of homogeneous equation 



(/ = 0), we have to obtain cq by other means, i.e. normalization of distribution. 

As for the name of the method, note that S n is expressed in terms of S n +i in the denominator, 
which can be in turn written in terms of S' n +2 and so on. It produces the continued-fraction structure 

s = m , P vtr - ( A - 5 ) 

11 t g2+ ... 

When the quantities in the recurrence relation are scalar, we call this the scalar continued frac- 
tion method. However, this method also applies to vectors recurrence relation, we then talk about 
matrix continued fraction. In such case, c n , f n and a n become vectors, while Q n and S n are matri- 



ces. The inversion in Eq. (A. 3) then becomes matrix inversion from the left (A/B = B 1 A). 



In fact, the recurrence relation can be treated as a set of linear equations, and solved by inverting 
a matrix of dimension N x N (for scalar recurrence relation). The operation involves complexity 
of 0(N 3 ). The use of continued fraction method reduces the complexity to O(N), and allows us to 
handle a much larger system of equations. 
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A. 2 Hubbard Operators 



Here we discuss some of the properties of the Hubbard operators X nm = \n){m\. They form a 
complete set, and one can think of X nm as a matrix with zeros everywhere, except 1 at the position 
(n, in). Some of the useful properties are 



Any operator can be expressed in terms of the Hubbard operators: 
A = } y A nm X nm , where A nm = (n\A\m). 



(A.6) 



Equal-time relation 



Commutator 



X n kXi m — OkiX nm . 



\X n }~^Xi m \ — ^k\X n m $mnXl]~. 



Adjoint 



{X nm )^ — X n 



Relation with the density matrix 

Xr[pA" nm ] — (X nm j — Pmn- 



(A.7) 



(A.8) 



(A.9) 



(A.10) 



Using the eigenstates of the system, the Heisenberg equation of the Hubbard operator becomes 

dX nm i 



dt h 



^nmX nm , 



(A.ll) 



where A nm = e n - e m . 



A. 3 Transition Rate 



To obtain the expression Eq. (5.34), one needs the identity 
x coth^Tx) 



f 

Jo 



(x 2 + A 2 )(x 2 + B 2 ) 



dx 



a^wH^)-^ 



7T 



1 



r AB(A + B) ' 



(A.12) 



in evaluation of the imaginary part. The Psi function, defined as tp{x) = ^lnT(x), has the following 
properties [22] : 



ip(x + l) 

m 

tp(x) 



i>(x) + 



i 



C where the Euler number C = 0.577215; 

t 



1 f°° 

2x~~ 2 J Jf + x 2 )(e 2 ^ - 1) 



dt. 



(A.13) 
(A.14) 

(A.15) 
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From the last expression with an integral, we obtain Eq. (A. 12) upon a partial fraction expansion 
of the denominator 



1 



1 



(x 2 + A 2 ){x 2 + B 2 ) B 2 -A 2 



1 



1 



x 



2 + A 2 x 2 + B 2 



(A.16) 



The real part and imaginary part of the transition rate are plotted in Figure |A.1| At low tem- 
perature, the real part decreases monotonically and is negligible at positive A, indicating that the 
process is dominated by de-excitation. At high temperature, the contribution from both excitation 
and de-excitation are important. 
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Figure A.l: In arbitrary unit, the transition rates at low temperature (top, (3 = 2) and high 
temperature (bottom, (3 = 0.05) with cut-off frequency Hojd = 10. 
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